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ABSTRACT 

Proteinaceous aggregation occurs through self-assembly— a process not entirely 
understood. In a recent article [1], an analytical theory for amyloid fibril growth 
via secondary rather than primary nucleation was presented. Remarkably, with 
only a single kinetic parameter, the authors were able to unify growth 
characteristics for a variety of experimental data. In essence, they seem to have 
uncovered the underlying allometric laws governing the evolution of filament 
elongation simply from two coupled non-linear ordinary differential equations 
(ODEs) stemming from a master equation. While this work adds significantly to 
our understanding of filament self-assembly, it required an approximate analytical 
solution representation. Here, we show that the same results are found by purely 
numerical means once a straightforward and reliable numerical solution to the set 
of ODEs has been established. 

Key Words: Master Equation, Proteinaceous aggregation, Self-assembly, 
Convergence Acceleration 

1. INTRODUCTION 

The prediction of proteinaceous aggregation may be the key to understanding the 
pathology of a host of degenerative transmittable diseases [2]. Whether a result of 
external initiation or a symptom of infection, protein misfolding seems to be a 
central element in prion protein disease progression. Fortunately, how 
macromolecules polymerize into long chains from monomer nucleation lends itself 
to simulation through kinetic equations that include self-assembly and disassembly 
[3]. As with any physical description and especially for biophysical systems, the 
specification of meaningful rate constants is essential for a successful prediction of 
chain length distributions. More importantly however, is the discovery of 
underlying allometric principles governing fibrillogenesis. This requires forming 
appropriate combinations of rate constants and investigating their invariance to 
self-assembly. That was the objective of a recent article in SCIENCE [1] 

[hereafter referred to as 1] where self- assembly promoting filamentous growth 
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was theoretically investigated in an attempt to establish fundamental scaling laws. 
In the article, based on a master equation characterizing nucleated polymerization 
through conventional kinetics, the authors remark 



"The lack of analytical solutions to such master equations has, 
however, represented a challenge to the quest to establish general 
principles and laws governing filamentous growth." 

The purpose of the present investigation is to demonstrate clearly that this is not 
the case and to indicate how the same fundamental scaling laws can be uncovered 
through purely numerical means. The necessary elements to do so are (1) an 
accurate numerical solution to a set of non-linear ODEs and (2) a straightforward 
dimensional analysis. While presented in the context of microscopic filamentous 
growth, on its own, the numerical approach represents a new way to consider a 
numerical solution to coupled non-linear ODEs. In particular, a highly accurate 
numerical solution algorithm is proposed rivaling the accuracy of any analytical 
solution, in particular that found in 1 . 

We begin with the master equation describing filamentous growth through 
monomer nucleation. As with all physical investigation, the formulation of a rate 
equation, balancing creation and destruction of the physical elements, is a central 

consideration. At this point, we refer to the presentation of 1 , where the following 
master equation for the filament length distribution function appears: 

&M=2kMt)f{t,j-i)- 

-2k t m(t)f{t,j)-k_{j-\)f(t,j)+ (1) 

oo 

+ 2k_Zf(tJ) + k n m(t) n S jn . 

i=j+\ 



The first term after the equals describes the creation of a filament of length j from 
the nucleation of monomers at either end of a filament of length y'-l. The second 
term represents the loss of a filament of length j as it grows to length y'+l through 
secondary nucleation. The third term specifies the possibility of a filament of 
length j breaking at any of its y'-l links. The fourth term refers to the contribution 
from all filaments of length greater than j breaking at either end to form a filament 
of length j, and the final term approximates the source of spontaneous and 
homogeneous monomer nucleation to length n c seeding the growth of all 
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subsequent filaments. is the monomer concentration found from the 

following moment equations: 



^ = k_ [m (t) + (ln c - l)P(t)] + k n m {tf + k_m tot 



(2a) 



where 



j= n c 

is the filament number density; and if 

00 

M(t)=Zjf(t,j), 

j= n c 

is mass concentration, then 

m(t) = m tot -M(t), (2b) 

where m tot is the total protein concentration. Equations (2a,b) are to be solved with 
initial conditions 

(2c) 

M(0) = M . 

The main thrust of our investigation is to obtain a highly accurate and efficient 
solution to Eqs (2) and in so doing enable a methodology for discovering 
underlying scaling laws governing filament growth. Our focus will first be on the 
development of an accurate numerical solution thus reducing numerical uncertainty 
and providing a benchmark for verification of more comprehensive numerical 
algorithms in the future. We emphasize accuracy through simplicity in order to 
make high accuracy results available to numerical experts and non-experts alike. 
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2. A FINITE DIFFERENCE NUMERICAL ALGORITHM 



A more convenient form of Eqs(2) emerges if one defines the moment vector 



>(<)' 

m(t) 



(3a) 



leading to 



where 



d ^ = A{y(,))y(t) + S(,), 



(3b) 



-k_(2n c -\) -k_ + k n m(t) 
n c (n c -\)k_ -2k + P(t)-n c k n m(tf 



n-\ 



(3c) 



and 



km tot 





The initial conditions then become 



(3d) 



*(0)- 



m, , - M,, 

tot . 



(3e) 



By uniformly partitioning the time interval into intervals of h = t j+l - t j and 
integrating Eq(3b) over an interval, we find 

yj« -yjA {y J+ i ) y J+ i + 4 ) ^ ] + \ V s ^ +s j1> ^ 



where a trapezoidal rule approximates the integrations. Hence, after matrix 
inversion, the solution is 
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y J+ i = 



I-^A J+1 (y J+1 ) 



However, because of the nonlinearity of the moments equation, Eq(4b) must be 
iteratively solved by lagging the matrix between times tj and t j+ \ 



The iterative solution therefore becomes 



y) + i 



i- h A J+1 

2 J 



(y'A) 





T k A 


]"{ 


I + -A. 




2 J 



m h 

y, +^ 



+ 



(5) 



To this point, the development has followed along classical lines from which we 
know the solution is of second to Eqs(4) is order [4] and will contain an iteration 
error as well. Our immediate goal therefore is to use both of these errors to our 
advantage to attain the highest possible accuracy. 

a. Romberg convergence acceleration 

A well known, but relatively little used method to obtain high accuracy for any 
finite difference solution is Romberg acceleration. The concept is the same as that 
behind Romberg integration but now applied to a sequence of approximations of a 
solution to coupled ODEs having known error variation. As already noted, the true 
solution takes the second order form 



^('y) = ^W + £ fl A.o A 

k=\ 



2k 



(6) 



where J/oW is the above finite difference approximation y^h) of Eqs(5) at 
iteration /, where j refers to the time tj of the initial grid. Then, eliminating the first 
term of the error series by considering y jfi {h) at the same time edit on a grid of half 
the original interval by evaluating y (h 1 2) , one finds 
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2 2 y. (h/2)-y. (h) 



2-1 



as the next highest order approximation (order 4). The true flux representation for 
the corresponding approximation becomes 



oo 



2 k 



k=2 



Continuing to eliminate higher orders sequentially by continually halving the grid 
gives the following recurrence relation for increasingly higher order 
approximations: 



yjA h ) s 



yjA h )=yj( h ) 

2^y hm _ 1 (h/2)-y j , m _ l (h) 



2 2m -1 



m = 1,2, 



(7) 



resulting in the convergence acceleration (or extrapolation) of the solution to zero 
discretization. The solution at the original edit, which is now t 2 „, . on the grid 

refined m-times, becomes 



oo 



(8) 



k-m+1 



It should be apparent that Romberg convergence acceleration applies only to the 
original time edit, which must be inherited by grid refinement. One way to ensure 
arbitrary edits is to perform grid refinement between these edits. 

Romberg acceleration requires that one generate the following sequence of finite 
difference approximations from Eqs(5) at each (inner) iteration /: 



y jfi (h/2 m ),m = 0,1,2,...., 



(9) 



to give the Romberg sequence for the recurrence of Eq(7). The Romberg 
acceleration, therefore, simply rearranges the original sequence into a more 
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efficiently converging one— hence the name convergence acceleration. To test for 
convergence therefore, we have the choice of convergence of the original sequence 



e - max 

k=\,2 



yjA hl2 ")-yiA hl2m ~ x ) 



< s 



(10a) 



or the Romberg sequence 



e„ = max 

R £=1,2 



y 'j,m,k V 1 ) y j,m-\,k 

yj, m ,k ( h ) 



(h) 



(10b) 



Note that we base convergence on the worst relative error of the two moments of y 
at each edit. 

In essence, we have redefined the meaning of a solution to be a sequence of 
solutions extrapolated to zero discretization error. No longer will a single 
discretization serve as the solution. 

b. Wynn-epsilon (We) convergence acceleration 

To treat the inner iteration resulting from the non-linearity, we use a Wynn-epsilon 
convergence algorithm [5] realizing that the iteration in Eqs(5) results in a 
sequence of solutions at each edit tj that presumably converges to the correct 
solution in the limit. In this way, we have created the sequence of approximations 

y l . during the iteration assumed to converge to the limit 



To accelerate convergence therefore, one can apply the Wynn-epsilon (We) 
convergence accelerator component-wise, which may result in a more rapidly 
converging sequence than the original. 

The We accelerator takes the following form 




(11) 
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sf=P] OTm) , l = 0,...,L 



-i -i 



, k = 0,...,L ; / = 0,...,L-£-l 



for either moment sequence Pj or m'j . The recurrence results in a tableau 



'0 



,(1) 

'0 



XL) 



,(0) 
'1 

,(1) 



,(0) 
'2 

,(1) 



,(i"2) 



,(0) 
'L-\ 

,(1) 



(12) 



where each element of an even column estimates the limit. Convergence comes 
from interrogation of last term of the even columns sf~'\ i = 0,2,...,2[Z/2] 



e We ~ 



£ {L~i) _ s (L-i-2) 



<s, i = 2,...,2[LI2\ 



(13) 



Actually, the We algorithm is a diagonal Pade approximant. 



3. NUMERICAL VERIFICATION 



The numerical solution of Eqs(3) is the application of We acceleration on the inner 
iterations of the nonlinear equation at each tj during the outer acceleration to a 
higher order accurate solution. In this section, we consider verification of the 
convergence acceleration algorithms thus far described as well as several 
additional convergence accelerations. Since the Romberg procedure requires a 
sequence of finite differences given by Eq(9), we can also attempt to accelerate 
this sequence through We acceleration. In addition, since the Romberg 
approximations are themselves a sequence, they can, in turn, be accelerated via the 
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We accelerator. Hence, the following four sequences are available for convergence 
interrogation: 



- the original 

- Romberg 



- We 



- We applied to Romberg. 



a. Verification by analytical solution 

If in Eq(3c), we let 




then 



-k -k +k 



n 







-k 



(14) 



n 



and Eqs(2) becomes linear. This case corresponds to no increase in filament length 
from nucleation. Filaments, however, are still able to break at their internal links 
as well as filaments of length larger than j can disassemble giving filaments of 
length j. However, only spontaneous growth from a single molecule is possible. 
For this case, the two differential equations for the moments are not only linear, 
but also they decouple into 



dP(t) 

dt 
dm{t) 



k_ \_m (t ) + P(t )] + kni (t ) + km 



k n m(t) 



tot 



(15a) 



dt 



with solution 



M(t) 
P{t). 




(15b) 
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In Figs. la,b, for the number and mass concentration of filaments for a total protein 
concentration m tot of 10/zM, we assume the parameters of Table 1. The approach 
of the moments to each other and to their asymptotic distributions is clear from 
Fig. la. A comparison of the convergence behavior in terms of the relative error 
with respect to the true solution for the four sequences is shown in the error chart 
of Fig. lb. The relative errors for each acceleration up to t = 15h, are presented 
with a time offset for clarity. Each horizontal line represents an increased 
discretization through seven levels of grid refinement. Most remarkably, as the 
original sequence reaches an accuracy of about 10" 5 in 7 grid refinements, the 
Romberg (R) and the We 



accelerated Romberg (We/R) achieve this in 3 levels; and for the We stand-alone, 
which overall was less impressive, in 4. This excellent performance clearly 
indicates the distinct advantage of acceleration. Indeed, we find near double 
precision machine accuracy. We next verify the algorithm including its nonlinear 
nature. 

b. Verification by Manufactured Solution 

A very powerful method of verifying the finite difference solution for nonlinear 
equations is through manufactured solutions [6]. In this method, we simply 
assume a solution, say P (V) and m (V), specify the source that gives this solution 

and finally run the numerical algorithm with the assumed source to see how well 
the assumed solution is reproduced. Here, we assume the solution to be the 

iterated analytical solution of 1 



Table 1 



Kinetic parameters for analytical solution 

k_=\o- 5 s~ l 



M =0.21 juM 
P =M /1380 M 
k n =22M_ M'V 1 
m M =10juM. 



P (t) = D + e Kt +D_e 



-Kt 




tot 



(16a) 
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Time Offset 

Fig. 1 . (a) Exact moment traces (b) True relative error chart 



M (t) = 



(16b) 
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and 



M {t) = m tot -m (t). (17d) 
Hence, the analytical moments solutions to Eqs(17a) are 

p(t)=p (t) 

m(t) = m (t). 

We are now in position to verify the full nonlinear algorithm. For this case, we 
assume the parameters of Table 2. Figure 2a shows the unbounded exponential 

Table 2 

Kinetic Parameters for the manufactured solution 
and converged numerical solution of Fig. 3 

k + =5x\0 4 M~V l 

£_=2xlO~V 

n c =2 

M =0juM 

P =0juM 

k n =2xl(TV 1 M- 1 

m tot = 5 /* M - 

behavior of the moments in dynamic equilibrium characteristic of homogenous 
nucleation. Again, the advantage of R and We/R accelerations are evident from the 
error chart of Fig, 2b. The standalone We acceleration however is again not as 
effective. For this case, we have used quadruple precision arithmetic to further 
highlight the advantage of convergence acceleration to 10" relative error. The 
inner iteration required a maximum of 10 iterations at each j, and we observe 
nearly the same error performance as found for the analytical benchmark. In 
addition, when the error chart, now for 8 grid refinements, is compared to that for a 
simple engineering estimate of the relative error between iterations, the error 
estimate always conservatively overestimates the true error. The computational 
effort required for both this and the first benchmark was less than 5 s on a Gateway 
1.4 GHz laptop. 
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With the algorithm now verified, we can apply it to the cases found in 1 to 
estimate the error of their "analytical solution" derived through fixed-point 
iteration. 



3. NUMERICAL CONFIRMATION 

In this section, we consider the results of 1 . The solution presented there, found 
by analytical iteration, 



P (0 = ~ e {lnrX)kJ Ei {-C + e Kt ) + B 2 e^- l)kJ (1 8a) 



2n -1 k 



M(t) 



m 



tot 



j n -1 



(18b) 



with 



c = ^P(0) ± M(0) ± ^-' 



K 



2m, , 2k 



'tot 



(18c) 



while expressed analytically, is by no means the true analytical solution. Actually, 
the analytical approximation (a more appropriate label) would eventually lead to 
the true analytical solution if the iteration were continued indefinitely. One of 
several obvious flaws in this representation is that the solution for polymer mass 
concentration M does not, in general, satisfy the initial conditions since for t = 



M(0) = 



m 



tot 



M(0) ' 

1m, , 
-e "* 



Only for a vanishingly small initial condition will this relation be satisfied. 
However, it is certainly an adequate approximation for an initial monomer 
concentration near m t , . 

tot 

Figure 3 shows the trace of the polymer concentration for the kinetic parameters of 
Table 2 for both the converged numerical solution and the analytical 

approximation of 1 . From this figure (also in 1 compared to a different numerical 
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Time (h) 




algorithm), the two solutions do not coincide on the steep rise of the sigmoid 
distribution. A closer examination however, reveals, in addition, the asymptotic 
values differ by a relatively small amount— 4.999997 compared to 5.000000 for 
the analytical approximation. While this is physically insignificant, it is not so 
when comparing two supposedly highly accurate solutions, which, in principle, 
should agree asymptotically to at least 7 digits. What is causing this difference 
then? 

To find the asymptotic distribution one sets the derivatives in Eqs(2) to zero giving 

= -k_ [m„ + {2n c ] + k n m n < + k_m tot 
= -2[k + m x -n c (n c -\)k_l2\P„ -njc n < 

where P m and m m are the asymptotic values if they exist. It is a simple matter to 
show 

p = m tOt ^ , K < 



2n c - 1 2n c - 1 k_ 2n c - 1 
2knK< +1 ~ [ 2 Kk- ~ n]k n k_ ] ml - (20) 
- \n in -\)k 2 +2k M k m, t \m^ - 2k 2 m, . = 

|_ c V c ) - + - tot J CO -tot 

For n c = 2, the later equation reduces to a cubic polynomial whose solution is 
known explicitly. In any case, assuming that « m tot , the second of Eqs(20) 
gives 

m ao U t 1 □ " (21) 

M + K - 



yielding 



M ro □ 6//M, 



which agrees exactly with the numerically converged solution at t = 60h and 
beyond. This simple exercise provides further confidence in the convergence of 
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the proposed algorithm and highlights the inaccuracy of the analytical 
approximation. 

One issue not addressed in 1 is the error in the lag phase r lag resulting from the 

difference of the analytical approximation and numerical solution at the inflection 
of the sigmoid response curve. r lag is a measure of the time at onset of rapid 

filament growth and is the time intercepted by the tangent line through the 
inflection point. Since the onset of rapid growth depends on this slope, the 
difference between the analytical approximation and the true solution should play a 
significant role in any analysis. An example of the inflection point, or maximum 
of the first derivative of the monomer population, is shown in Fig. 3b. From Fig. 
4, where the phase lag for 2401 random variations of the kinetic parameters each 
over two decades is shown, it seems that r lag is overestimated by the analytical 

approximation (as expected because of its steeper slope). Fig.4 required an 
execution time of about 25 min for a relative error of 10" 5 . 
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Numerically converged solution 
Analytical approximation 



The variation between the two solutions is also apparent in comparisons with 
experimental measurements. The polymer concentrations shown in Fig. 5 
represent such a comparison with the kinetic parameters of the analytical 
approximation chosen to best represent the data. In general, the numerical solution 

is mostly within the errors bars of the data (see 1 ) and in some instances provides a 
more appropriate fit. The parameters for these cases are given in Tables 3a,b. 

Table 3a 
Parameters for Figure 5a 

k + =2.9xl0 4 M-V 
^_ =2.1xl0" 9 s _l 
n c =2 

M(0) = 0.21//M 
k n /k_=22.SM~ l 

Table 3b 
Parameters for Figure 5b 

k + =2.9xl0 4 M-V 1 
k_ = 2.1xl(TV 

n c =2 

M(0)/P(0) = 1380 
k n /k_=22.SM~ l 
m tot = 98//M. 

4. DIMENSIONAL ANALYSIS 

With the numerical algorithm verified and the numerical results of 1 confirmed, 
we are ready to use it in search of the scaling laws governing amyloid filament 
growth. 

As already noted, the authors of 1 claim that to establish allometric scaling laws 
for proteinaceous filament growth, analytical representations (or approximations as 
they have developed) are required. Hence, the claim is that numerical solutions by 
themselves are limiting. This is certainly true if we do nothing more than just find 
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an accurate finite difference scheme as we have done. As we shall show however, 
coupling a confirmed numerical solution to a dimensional analysis can potentially 
reveal the desired scaling laws. 



With knowledge only of modeling parameters, independent and dependent 
variables, the solution to Eqs(2) (for n c = 2) can be expressed as 

y(t) = F[k_,k + ,k n ,m tot ;t,P(t),M(t)~] (22) 

If we examine the units of all quantities, one finds they are combinations of only 
two fundamental units— moles (M) and time (T). In terms of these units, the listing 
for the parameters and variables in Eq(22) are given in Table 4. According to the 
Buckingham Pi theorem [7], since there are 7 parameters and variables and two 
fundamental units, we can express the solution in terms of five dimensionless 

Table 4 

Dimensions of parameters and variables in Eq(22) 



IK] 


=T -i 


[*-] 


= T A M 


[K] 


= T l M 


[ m tot] 


= M 


w 


= T 




= M 


[M(t)} 


= M 



quantities in the form 

^ = G[x l J = \,...,5] (23) 

m ,ot 

where n i is a dimensionless variable. To obtain a convenient set of these 

variables, we begin by specifying the following cumulative dimensionless variable 
containing all possible choices: 
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7r(a,b,a,j3,y;7r p ,7r M ) = 7r(a,b,a,j3,y) 7r p (t) x M (t), (24a) 

where 

x(a,b, a, f3, y) = k a k b + K<t r ■ (24b) 
Note that we anticipated two of the dimensionless variables as 

, n. Pit) , . Mit) 

> —■ ( 24c ) 

m tot m tot 

We must now determine the remaining three groups from the dimensional 
equivalence 

[x(a,b,a,P,y)\ = T r<a+b+a) M p - b - a . (25) 
Therefore, for 7r(a,b,a,j3,y) to be a dimensionless quantity 

y-(a + b + a) = 
P-b-a=0 

or expressing two of the five unknowns in terms of the remaining three gives 

a = y-P 
b = p-a; 

and Eq(24b) becomes 

x(a,b,a,j3,y) ^x(a,j3,y) = k^k^K^tf ■ (26) 

Since the specification of the dimensionless groups in not unique, it is convenient 
to specify all possibilities of a given class and choose an appropriate set. For this 
reason, we consider all dimensionless variables to powers of 1 or - 1 by letting a, ft, 
and /cycle through and 1. This produces two sets of results as shown in Table 5. 
The choice of which three of the six groups is primarily one of convenience. 
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If additional values of a, /?, and y are assigned, additional groups are formed as 
products of the groups already found. In particular, for the combination (0,1,2), we 
find 

n n (t) = k k + m t J 2 = n\n 1 = kI2, (27) 

which plays an important role in the analysis of 1 and is called the rate of 
multiplication of the filament population. 

Table 5 



I 


a,P,7 




1 


0,0,1 


kj 


2 


0,1,0 


K 

tot k_ 


3 


1,0,0 


K 
K 


4 


0,1,1 


k + t 


5 


1,1,0 


tot j 

k_ 


6 


1,1,1 


kjn t J 

n tot 



The truth of the dimensional analysis is immediately obvious from the analytical 
approximation of 1 since, in terms of the dimensionless groups, Eqs(18) become 

n p (t) = -- e^Ei^-C^^ ) + B 2 e-^ {t) 

3 V 2;r 2 ' (28a) 

7T M (0 = 1- zM-C + e M ^ + C/ ff ' W ^ + 7T 2 7T 3 ) 

with 
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C ± (n p (0),7T M (0),7T 2 ,7T 3 ) = 71 p (0)^2^±-7T M (0) ±-7U 2 7C z 

" (28b) 

B 2 (7r p (0),7r M (0),7r 2 ,7r 3 ) = 7r P (0)-- + n =Ei(-C + ). 

With the dimensionless groups as our guide, we now search for invariances in 
experimental data. For example, say we have the data shown in Fig. 6 [8]. We 

/ \l/2 

first consider 7T 5 and ti 1 ( r lag ) since these groupings contain all four parameters 

and time. Figure 7a is a plot containing the 2401 random values of the kinetic 
parameters used to generate Fig. 4 and these two groups. Remarkably, as predicted 
in 1 , we find an almost perfect linear correlation 

\nx 5 = -1.132-1.414^^") 



1.2 




5 10 15 20 25 30 

Time (h) 

Fig. 6 Approximations to Ferguson data [8]. 
which when manipulated gives 
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T lag 



yjlk + k_m 



:[ln(l/C + )- 1.825] 



(29) 



tot 



and is the scaling found in 1 only differing in the subtracted constant (1.825 vs. 
1.718) to give a longer lag time. Since this expression comes from an accurate 
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Fig. 7a. Scaling law for lag phase. 

numerical solution, Eq(29) represents a more accurate scaling law and clearly 
demonstrates that a purely numerical solution can indeed provide such. Finally, 
when the experimental data, numerically modeled in Fig 6, is added in addition to 
the data for Fig. 3, the correlation is confirmed. 

There is additional information to be found from our numerical investigation 
however. The normalized maximal growth v m (defined as the derivative at the 
inflection point normalized to the total protein content) when multiplied by the lag 

phase also forms a dimensionless variable. When plotted against 7£ 7 [T lag j as 
shown in Fig. 7b, again a nearly perfect correlation results giving the scaling law 
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„ /—— 5.548xl0" 4 

v m = 0.3 1 ^2k + km tot + , (30) 

which is again confirmed by experimental data. Hence, once the lag phase is 
found from Eq(29), we also find the maximal growth from Eq(30). Finally, when 
Eq(29) is introduced into Eq(3), \ m is observed to be proportional to k, which is 

consistent with 1 . 




5. FINAL COMMENTS 

We have developed a new numerical solution for the moments equations of 
amyloid fibril growth resulting from a master equation. The solution features high 
accuracy through convergence acceleration using both Romberg and Wynn-epsilon 
procedures to extrapolate a simple finite difference approximation to nearly zero 
discretization error. The numerical strategy is verified by comparison to an 
analytical solution of a related linear set of ODEs and to a manufactured solution 

based on the first iteration of the analytical approximation found in 1 . While the 
numerical algorithm in itself represents a new approach to solving ODEs, its value 
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is further enhanced when coupled to a dimensional analysis. In this way, a 
dimensional analysis provides a guide to the underlying evolution of fibril growth 
from which the fundamental scaling laws emerge. Indeed, we confirm the scaling 

laws found in 1 to a high level of confidence since they result from a more 
accurate numerical solution to the governing moments equations. As a final 
example of this, another scaling law emerges when comparing dimensionless 
groups 7t\ and x 2 as shown in Fig. 8. This scaling is not as strong as the two above 
however. Also indicated are the experimental data of Fig, 6 and the analytical data 
of Fig. 3. While not in complete agreement with the correlation, the law is 
indicative of the trend. 

The significance of this work is the demonstration that an analytical solution is not 
required to discover scaling laws. All that is needed is a reliable numerical 
solution and a dimensional analysis. Most likely this will be the preferred 
(probably the only) approach when reaction rates are dependent on the component 
concentrations which seems to be the direction investigations are headed. In 
addition, it will be the approach to the solution of the master equation itself, which 
is now relatively easily within reach using convergence acceleration. 



1e-2 




1e-6 -I 1 1 1 1 1 

1e+6 1e+7 1e+8 1e+9 1e+10 1e+11 1e+12 

Fig. 8. An additional scaling law uncovered by the numerical 
formulation 
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